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Abstract — In this paper we investigate the continuum limits 
of a class of Markov chains. The investigation of such limits 
is motivated by the desire to model very large networks. We 
show that under some conditions, a sequence of Markov chains 
converges in some sense to the solution of a partial differential 
equation. Based on such convergence we approximate Markov 
chains modeling networks with a large number of components 
by partial differential equations. While traditional Monte Carlo 
simulation for very large networks is practically infeasible, partial 
differential equations can be solved with reasonable computa- 
tional overhead using well-established mathematical tools. 

Index Terms — Continuum modeling, Markov chain, partial 
differential equation, large network modeling, wireless sensor 
network. 



I. Introduction 

NETWORK modeling is an important tool in the analysis 
and design of networks. Many network characteristics 
of interest can be modeled by Markov chains, where Monte 
Carlo simulation has been the traditional approach [1]. With 
the enormous growth in the size and complexity of today's 
networks, their simulation becomes more computationally 
expensive in both time and hardware. Some effort has been 
made to exploit the computing powers of distributed computer 
networks, such as parallel simulation techniques, where the 
number of processors needed in the simulation increases with 
the number of nodes in the network 13, Q. However, for 
networks involving a very large number of nodes, Monte Carlo 
simulation eventually becomes practically infeasible. 

In this paper we address this problem by focusing on the 
global characteristics of an entire network rather than those 
of its individual components. The idea is to approximate 
the underlying Markov chain modeling a certain network 
characteristic by a partial differential equation (PDE). 

As a concrete familiar example, which we present in Sec- 
tion HIl consider multiple i.i.d. (independent and identically 
distributed) random walks of M particles on a network con- 
sisting of TV points. For any vector x, let x T be its transpose. 
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Let the Markov chain modeling the network characteristic 
be X N (k) = [XAr(fc,l),...,XAr(>,iV)] T € R N , where 
Xjy(k,n) is the number of particles at point n at time k. If 
we treat N and M as indices that grow, this defines a family 
of Markov chains indexed by N and M. We show that as 
M — > oo and TV — » oo, Xjv(fc) converges in some sense to its 
continuum limit, a deterministic function with continuous time 
and space variables. Under certain conditions, it is possible to 
characterize such a function as the solution of a PDE [4]-|6l. 
This itself is not a new result, but helps to illustrate our aim. 

Indeed, our development here is motivated by the network 
modeling strategy in [7] and the need for a rigorous description 
of its underlying limiting process. We illustrate in Section iHll 
the convergence of the sequence of Markov chains to the PDE 
in a two-step procedure. Suppose the evolution of Xj^(k) is 
governed by a certain stochastic difference equation with a 
"normalizing" parameter M, Let X]y(k) be the normalized 
deterministic sequence governed by the corresponding "ex- 
pected" and deterministic difference equation. First, we show 
in Section IIII-BI that Xfj(k)/M is close to ccjv(fe), in the 
sense that as M — > oo, both their continuous-time extensions 
converge to the solution of an ordinary differential equation 
(ODE). Second, we show in Section IIII-CI that as N — > oo, 
xjv(fc) converges to the solution of a PDE. Therefore, as 
M -> oo and N — > oo, X N (k)/M converges to the PDE 
solution. 

Our procedure provides an approach to approximating 
Markov chains that model large networks by PDEs. PDEs are 
widely used to formulate time-space phenomena in physics, 
chemistry, ecology, and economics (e.g., ll8ll- lflTll ). and there 
are well-established mathematical tools for solving them such 
as Matlab and Comsol, which use finite element method lfl2ll 
or finite difference method fl3l . In contrast to Monte Carlo 
simulation, our approach enables us to use these tools to 
greatly reduce computation time, which makes it possible 
to carry out the analysis, design, and optimization for very 
large networks. We present in Section [IV] an example of 
the application of our approach to the modeling of a large 
wireless sensor network. In this example, we derive an explicit 
nonlinear diffusion-convection PDE, whose solution captures 
the dynamic behavior of the data message queues in the 
network. We show that although the PDE approximation takes 
only a tiny fraction of the computation time of the Monte 
Carlo simulation, there is a strong agreement between their 
simulation results. 

Continuum modeling has been well-established in fields 
such as physics, mechanics, transportation, and biology (e.g., 
(Si-El)- Its applications in communication networks, how- 
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ever, are relatively new and rare. Among these, to our best 
knowledge, our approach is the first to address the time- 
space characteristics of communication networks with a large 
number of nodes. In contrast, for example, [18|-[20| deal with 
networks with heavy traffic instead of large number of nodes; 
lED . [22] present scaling laws of the network traffic without 
characterizing the actual traffic over time and space; and |23|, 
[24 1, which use mean field methods, only keep track of the 
statistical features of the networks such as the fraction of nodes 
in each network state. 



II. Continuum Limit of Multiple Random Walks 

In this section we present an illustrative example of approxi- 
mating multiple i.i.d. random walks by a PDE. First consider a 
single random walk on a one-dimensional network consisting 
of N points uniformly placed over V = [0,1], as shown in 
Fig. [U Hence the distance between two neighboring point is 
ds = 1/(N + 1). At each time instant, the particle at point 
n, where n = 1 , . . . , N, randomly chooses to move to its 
left or right neighboring point with probability P r (n) and 
Pi(n), respectively. Let the length between two time instants 
be dt = 1/M. We set dt = ds 2 , which is a standard time-space 
scaling approach to ensuring the convergence of the difference 
equation to a PDE. We assume a "sink" boundary condition, 
i.e., the particle vanishes when it reaches the ends of T> (though 
"walls" at the boundary are equally treatable). 
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1 n-l n n+1 N N+l 

Fig. 1. An illustration of a one-dimensional single random walk. 

Now consider M random walks on the same network, 
where the particle in each random walk behaves independently 
identically as in the single random walk described above. 
Let Bi(k,n) be the Bernoulli random variable representing 
the presence of the ith particle at point n at time instant 
k, where k = 0, 1, . . . and n = 1, . . . , N. Define Bi(k) = 



[Bi(k,l),...,Bi(k,N)] J 



piV 



According to the behavior 



of the particle in the single random walk, for i — 1, . . . , M, 

Bi(k + l,n) - Bi(k,n) 

Bi(k,n — 1), with probability P r (n — 1); 

Bi(k, n + 1), with probability Pi(n + 1); 

—Bi(k,n), with probability P r (n) + P/(n); 

0, otherwise, 

where Bi(k, n) with n^Oorn^A^ + l are defined to be 
zero. Let the function Fn(x, U(k)), where U{k) are i.i.d. and 
do not depend on x, be such that for i — 1, . . . , M, 

Bi(k + 1) = Bi(k) + F N (B t {k), U{k)). (1) 

Then for x = [x±, . . . , xn] t , the nth component of 



Fn(x, U(k)), where n = 1, . . . , N, is 

x„_ 1; with probability P r (n — 1) 

x n+ i, with probability Pi(n + 1) 

—x n , with probability P r (n) + Pi(n) 

0, otherwise 



(2) 



where x n with n^Oorn^A^ + l are defined to be zero. 

Let Xpj(k, n) be the number of particles at point n at time 
k. Then 



M 



X N (k, n) = Bi(k, n). 



(3) 



Define X N (k) = [X N (k, 1), . . . , X N (k, N)] T , which forms a 
discrete-time Markov chain with state space R N . Since _FV is 
linear, it follows from (f3]) that 



X N {k + 1) = X N (k) + F N (X N (k), U(k)). 



Let 



f N {x) = EF N {x,U{k)), x€ 



i>/Y 



It follows from (|2]i that for x — [xi, . . . , xn] T , the nth 
component of Jn(x), where n = 1, . . . , N, is 

P r (n - l)x n -i + P{n + l)x n+1 - (P r (n) + P(n))x n , (4) 

where x n with n^Oorn^A^ + lare defined to be zero. 
By (Q]l and the linearity of Fjy, for i = 1, . . . , M, 



EB t {k + 1) = EBiik) + f N (EBi(k)). 



(5) 



Notice that, since the random walks are i.i.d., EB L does not 
depend on i. Define a deterministic sequence XN(k) by 



x N (k + 1) = x N (k) + fN(x N (k)), 



(6) 



where 



xn(0) = — Tj ' a ' s ' ( a l most surely). (7) 

We seek to approximate Xjv(fc) by a continuum model, 
where the time and space indices k and n are made continuous 
as N — > 00 and M — > 00 in the following two steps: First, 
define 

X N ([Mt\) 



X Q N(t) 



M 



the continuous-time extension of X^(k) by piecewise- 
constant time extensions with interval length dt = 1/M and 
scaled by 1/M. Second, define X p ff(t, s) to be the continuous- 
space extension of X n {t) by piecewise-constant space exten- 
sions on V with interval length ds. Notice that as N — > 00, 
ds — > 0. Thus X p n is the continuous-time-space extension of 
Xjs[{k). Similarly, define x N{t) — XN([Mt\), the piecewise- 
constant continuous-time extension of xpf(k), and x P N(t,s), 
the piecewise-constant continuous-space extension of x ^(t). 
Thus x p m is the continuous-time-space extension of xjv(fc). 

Now we show that for M sufficiently large, X p n, the 
continuous-time-space extension of Xpj(k), is close to x p n, 
the continuous-time-space extension of xjv(fc). By OJ and the 
strong law of large numbers (SLLN), for each k, 



lim 



M- 



X N (k) 
M 



= EBi(k) a.s. 
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By this and 0, 



lim x N (0) = EBi(0) a.s. 



By (0 and (0, xpj(k) and EBi(k) satisfy the same difference 
equation. Then we have for each k, 



lim XM{k) — EBAk) a.s. 

M— too 



Hence for each k, 

lim = x N (k) ,s. 

Therefore, AT jv and x n are close for large M in the sense 
that 

lim \\X oN {t) - x oN {t)\\W =0a.s, (8) 

M^oo 



where 



|(JV) 



is the oo-norm on 



pJV 



Note that 



ll-Xpiv(-»*) - afpJvC-**)!!^ = \\X oN - x oN \\£\ 

2/ is the oo-norm on MP, the space of functions 
of T> — > R. Then by (0, X p jv and a; p Ar are close to each other 
for large M in the sense that 



where 



M 



lim \\X pN (;t) - x pN (;t)\\W =0 a.s. 



(9) 



Therefore, we can approximate X p n by a; p Ar for M suffi- 
ciently large. 

Next we show that as N — > oo, aip^r satisfies a certain PDE 
that is easily solvable. By (0 we have for n — 1, . . . , N, 

x^(k + 1, n) — XN{k,n) 

= P r {n - l)x N (k,n- 1) + + l)xjv(fc, n + 1) 
- (P r (n) + Pi{n))x N (k,n), 

where xjv(fc, n) with n < or n > N + 1 are defined to be 
zero. Assume Pi(n) = pi(nds) and P r (n) = p r (nds), where 
Pi(s) and p r (s) are real-valued functions defined on V. Then 
by the definition of x v n, it follows that for s £ P and t > 0, 

£ p at(£ + s ) _ %pN(t, s) 

= p r (s — ds)x p N(t, s — ds) + pi(s + ds)x p N(t, s + ds) 

- (Pr(s) +Pl{s))x pN (t, s). (10) 

To ensure a finite non-degenerate limit, we assume 

Pi(s) = b(s) + ci(s)ds and p r (s) = b(s) + c r (s)ds. 

Define c = q — c r . We call b the diffusion coefficient and c 
the convection coefficient, for a greater b means more rapid 
diffusion and a greater c means a larger directional bias. 
Assume that b 6 C 2 and c 6 C 1 . Assume that x p n is 
twice continuously differentiable in s. Put into ( TTOb the Taylor 
expansions 



X p N{t, s ± (is) = X p N(t, s) ± 
<9 2 a;pjv 



<9;r 



piV 



(t, s)ds 



+ 



ds 2 



8s 

(t,s)— + o(ds 2 ), (11) 



ds 2 



6(s ± ds) = b(s) ± b s (s)ds + b ss (s)^- + o(ds 2 ), (12) 



and 



c(s ± (is) = c(s) ± c s (s)ds + o(ds), 



(13) 



where a single subscript s represents first derivative and a 
double subscript ss represents second derivative. Then we have 



x P N(t + dt, s) — x p iv(t, s) = b(s)- 



(2b s (s) + c(s)) 



Ox 



ds 2 



(t, s)ds 2 



pN 



ds 



(t, s)ds 2 



+ Q> ss (s) + c s {s))x pN (t, s)ds 2 + o(ds 2 ). 
Divide both sides of (fT~4T > by dt = ds 2 and get 

x P N (t + dt,s) — x pN (t,s) 



(14) 



6(s)- 



dt 

2 X p N 



ds 2 

+ (b ss (s) - 



(t,s) + (2b s (s) + c(s)) dXpN 



ds 



(*,*) 



,(s))x pN (t, i 



o(ds 2 



ds 2 ' 

As N — » cxd, ds — > 0, and hence dt = ds 2 0. Assume 
that x p N is continuously differentiable in t. Then by taking 
the limit as N — > oo and rearranging, we get a PDE that x p n 
satisfies: 

x pN {t,s) = ^ (b(s)^^-(t,s) 



ds 



d 



+ -g^((b s (s) + c{s))x pN (t, s)), 

for t > and seP, with boundary condition x P N{t, s) = 0. 

As — s> oo, dt = ds 2 — > 0, and hence M = 1/dt = 
1/ds 2 —> oo. Then by (0, for N sufficiently large, X p n, 
the continuous-time-space extension of ATjv(fc), is close to 
x P N, the continuous-time-space extension of XN{k). There- 
fore, we can approximate XN(k) by the solution of the 
above PDE called the one-dimensional diffusion-convection 
equation, which can be easily solved l25l . Note that our 
derivation here differs from that of the well-studied Fokker- 
Planck equation (also known as the Kolmogorov forward 
equation) [26], whereas the latter originates from the study 
of the probability density of a Wiener process. 

This motivational example raises some questions that must 
be answered by the convergence analysis of the underlying 
limiting process. First, general networks may exhibit more 
complex behaviors. For example, Fn might no longer be lin- 
ear; and SLLN might not apply in many scenarios since node 
behaviors are not necessarily i.i.d. Specifically, the analysis 
above does not apply to the network Markov chain in Q. 
To find the conditions under which (0 holds in more general 
setting, in Section lTlI-Bl we apply Kushner's weak convergence 
theorem in [4] to a more general class of systems modeled by 
Markov chains. Moreover, we need to show in what sense and 
under what conditions X p n converges to the solution of the 
PDE. We analyze such convergence and provide its sufficient 
conditions in Section UlI-CI 

III. Continuum Limits of Markov Chains 

In this section we analyze the convergence of a sequence 
of Markov chains to the solution of a PDE in a two-step pro- 
cedure. We provide sufficient conditions for this convergence. 
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A. General Setting 

Consider N points placed over a Euclidean domain D 
representing a spatial region. We assume that these points form 
a uniform grid, though our approach can later be generalized to 
nonuniform cases. We will refer to these N points in V as grid 
points and denote the distance between any two neighboring 
grid points by ds^. 

Consider a discrete-time Markov chain 

X N (k) = [X N (k,l),...,X N (k,N)] T (15) 

with state space R w . Here Xisr(k,ri) is the real-valued state 
of point n at time k, where n = 1, . . . , N is a spatial index 
and k = 0, 1, . . . is a temporal index. 

Suppose that the evolution of X^(k) is described by the 
stochastic difference equation 

X N (k + l)=X N (k) + F N (X N (k)/M,U(k)), (16) 

where U(k) are i.i.d. and do not depend on the state .XV(fc), 
M is a "normalizing" parameter, and Fjy is a given function. 
Let 

f N (x) = EF N (x,U(k)), x£R N . (17) 
Define a deterministic sequence xjv(fc) by 

x N (k + l)=x N (k) + -^f N {x N (k)), (18) 

where xjv(O) = Xjv(0)/M a.s. In the next subsection, we 
show that under certain conditions, X^(k)/M and XAr(fc) are 
close in some sense. 

B. Convergence to ODE 

Let X N(t) be the continuous-time extension of Xw(k) by 
piecewise-constant time extensions with interval length 1/M 
and scaled by 1/M, i.e., for arbitrary t£R, 

X oN (t) = X N ([Mt\)/M. (19) 

It follows that for each k, X oN (k/M) = X N (k) /M. Similarly 
we define x jv(i), the continuous-time extension of xjv(fc) by 

x oN (t)=x N ([Mt\). (20) 

For fixed f N > 0, let D N [0,f N ] be the space of Re- 
valued Cadlag functions on [0,7V], i.e., functions that are 
right-continuous at each t £ [0,2V) and have left-hand limits 
at each t £ (0,2V]- As defined in ( [T9j > and < f2Qb respectively, 
both X oN (t) and x oN (l) with I £ [0,2V] are in D N [0,f N ]. 
Since both X iy(t) and x jv(i) depend on M, each one of 
them forms a sequence of functions in D [0, TV] indexed by 
M = l,2,.... _ 

Define the oo-norm || • ||^ on D N [0, Tjy], i.e., for a; G 
^[O.TV], 

||x[[W = max sup |x n (t)|, 
n=i,...,N te[0i f N] 

where x n is the nth components of x. A sequence of functions 
xm £ D N [0,Tjf] is said to converge uniformly to a function 
x £ D N [0,f N ] if as M -> oo, \\x M ~ x\\ ( £ -> 0. In this 

paper, we use the notation "=>" for weak convergence and 

p 

" — for convergence in probability. 



Let fff be defined as in (fTTT i. Now we present a lemma stat- 
ing that under some conditions, as M — > oo, X n converges 
uniformly to a limiting function y, the solution of the ODE 
V = /jv(y)> on [0,2V], and X at converges uniformly to the 
same solution on [0,7V]- 

Lemma 1: Assume: 
(la) There exists an identically distributed sequence {A(fc)} 

of integrable random variables such that for each k and 

x, \F N (x,U(k))\ < A(fc) a.s.; 
(lb) the function iV(x, U(k)) is continuous in x a.s.; and 
(lc) the ODE y = fpfiv) has a unique solution on [0,7V] 

for any initial condition y(0). 

Suppose that as M — > oo, 

X oN (0) A y(0) and 

Then, as M —> oo, 

\\XoN ~ y\\£ A and \\x oN - y\\ W -> 

on [0,2V], where y is the unique solution of y = /iv(j/) with 
initial condition y(0). 

To prove Lemma Q] we first present a lemma on weak 
convergence due to Kushner J4). 

Lemma 2: Assume: 
(2a) The set 

{\F N (x,U{k))\ : k > 0} 

is uniformly integrable; 
(2b) for each k and each bounded random variable X, 

lim E sup \F N (X, U{k)) - F N {X + Y, U(k))\ = 0; 

\Y\<5 

and 

(2c) there is a function /V(-) [continuous by (b)] such that 

as tl y oo, 

1 " 

-J2F N (x,U(k)) ^ f N (x). 

k=0 

Suppose that y = /jv(y) has a unique solution on [0,7V] 
for each initial condition, and that X o jy(0) 2/(0). Then as 

M oo, 

[[^oiV-2/ll^ ^0on [0,2V]. 

We note that in Kushner' s work, the convergence of X n to y 
is stated in terms of Skorokhod norm [4|, but it is equivalent 
to the oo-norm in our case where the functions are defined on 
finite time intervals l27l . 

We now prove Lemma Q] by showing that the assumptions 
(2a)-(2c) in Lemma [2] hold under the assumptions (la)-(lc) 
in Lemma Q] 

Proof of Lemma Q} 

1) Since X(k) is integrable, as a — > oo, 

^|A(fc)|l{|A(fc)|>a} 0, 

where 1a is the indicator function of set A. By Assump- 
tion (la), for each k and x, 

F N (x,U{k)) < A(fc) a.s. 
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Therefore for each x and a > 0, 

E\F N (x,U(k))\l {\F N (x,U(k))\>a} 

< £;|A(fc)|l{| Fjv ( Xii 7(fe))| >a } 

< E\\(k)\l {lml>a} . 
Hence as a — > oo, 

$upE\F N (x,U(k))\ln FN , XtU , m>a \ -> 0, 

fc>0 

i.e., the family { |iV (a;, [/(£;)) | : k > 0} is uniformly 
integrable and Assumption (2a) holds. 
2) By Assumption (lb), Fpf(x,U(k)) is continuous in x 
a.s. Then for each bounded X and each k, 

lim sup \F N {X, U{k)) ~F N {X + Y, U(k))\ = a.s. 



<s->o 



y|<<5 



By Assumption (la), for each x and each fc, there 
exists an integrable random variable A(fc) such that 
|Fjv(x, ?/(&)) | < A(fc) a.s. It follows that for each 
bounded X, each k, and each Y such that \Y\ < S, 

\F N {X, U(k)) - F N (X + Y, U(k))\ 

< \F N {X, U(k))\ + \F N (X + Y, U[k))\ < 2A(fc). 

Hence for each 5, 



sup \F N (X,U(k)) - F N (X + Y,U(k))\ 

\Y\<8 



< 2A(fc), 



an integrable random variable. By the dominant conver- 
gence theorem, 

lim.E sup \F N {X,U{k)) - F N {X + Y,U{k))\ 

\Y\<8 

= E]xm sup \F N (X,U(k)) - F N (X + Y,U(k))\ 

\Y\<S 

= 0. 

Hence Assumption (2b) holds. 
3) Since U(k) are i.i.d., by the weak law of large numbers 
and the definition of Jn in ( TT7| >, as n — >• oo, 



1 " 



k=0 



In{x). 



Hence Assumption (2c) holds. 
Then, by Lemma [2] as M — > oo, \\X n 



y\\ 



(o) 



on 



[0, Tff], For each sequence of random processes {X n }, if A 
is a constant, X, 



M-¥oo, \\X oN -y\\$- 



A if and only if X n A. Therefore, as 
on [0,T/v]. The same argument 

oo, 
■ 



implies the deterministic convergence of x n- as M 
\\xon ~ -»• on [0,f N ]. 

Based on Lemma [T] we get the following lemma, which 
states that X n and x n are close with high probability when 
M is large. 

Lemma 3: Let the assumptions in Lemma [T] hold. Then for 
any sequence {Cn}, for each AT, and for M sufficiently large, 
we have 

P{\\X oN - x oN \\t ] > Gv} < l/N 2 on [0,f N \. 



Proof: By Lemma Q] for each N, as M —> oo, 
\\XoN - y\\£ A and \\x oN - -y on [0,2V]. 
By the triangle inequality 

ll^oAr - x off \\g < \\X oN - y\\g + \\x oN - y\t\ 

it follows that as M —> oo, \\X n — a; w||oo — > on [0,7V]. 
This finishes the proof. ■ 

Since X n and x n are the piecewise continuous-time ex- 
tensions of X^ and xn by constant interpolation, respectively, 
we have the following corollary. 

Corollary 1: Fix f N and let K N = [f N M\. Let the 
assumptions in Lemma [T] hold. Then for any sequence {Cat}, 
for each N, and for M sufficiently large, we have 



P 



max 

*=0 Kji 

n = l y -.- y N 



X N (k,n) 



M 



N 



< 



1 



We use Lemma [3] and Corollary Q] in the next subsection. 

C. Convergence to PDE 

In the last subsection, we stated conditions under which the 
continuous-time extensions of X^(k) and xj^(k) are close 
asymptotically (as M — > oo) with high probability. In this 
subsection, we further let N — > oo and state conditions under 
which xjv(fc) is close asymptotically to the solution of a 
PDE. This leads to the convergence of X^(k)/M to the PDE 
solution as M — > oo and N — >• oo. 

Assume that the domain V introduced in Section IIII-AI 
is compact and convex, and let u; : T> — > 1 be in C 2 . 
Given a fixed N, let Vjsr be the set of the N grid points 
in T>. Let y^ be the vector in composed of the values 
of w at the grid points Ujv(n) 6 VV, n = 1, ...,N, i.e., 
y N = [w(v N (l)), w{v N (N))] T . 

Given s £ T>, let {sn} CPbea sequence of grid points 
in V such that as N ~ > oo, sn — > s, where for each iV, 
sjv is a grid point in Vn- Let /at(j/at, s^v) be the component 
of the vector /jv(2/jv) corresponding to the location sjv- For 
example, for = 5, if S5 = Ws(4) in V5, then /s(j/5,S5) is 
the 4th component of the vector /s(j/5). 

Assume that there exist sequences |(5at}, {/3jv}, {7Af}, and 
{pn}, functions / and h, and < c < 00, such that as N — > 
00, <5 at -> 0, 5n/ Pn — > 0, 7jv — ^ 0, pjv — ► 0, and: 

• for any sn such that sn — >• s, where s is in the interior 
of T>, there exists a sequence of functions cf>N ■ — > M. 
such that 

fN(yN,s N )/$N = f{s N ,w(s N ),Vw(s N ),V 2 w(s N )) 
+ <Pn(sn), (21) 
and for sufficiently large, 

|0iv(sjv)| < C7 W ; (22) 

and 

• for any s n such that sn — >• s, where s is on the boundary 
of 2?, there exists a sequence of functions tp N ; V — > R 
such that 

fN(yN,s N )/ '(3 N = h(s N ,w(s N ),Vw(s N ),V 2 w(s N )) 
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<Pn(sn), 



(23) and 



and for TV sufficiently large, \(Pn(sn)\ < cpjv. 
Here, V l w represents all the ith order derivatives of w, where 
i = 1,2. 

These assumptions are technical conditions on the asymp- 
totic behavior of the sequence of functions /jy. The basic idea 
is that /wiUNiS^) is asymptotically close to some function 
of terms that look like the right-hand side of a time-dependent 
PDE. Typically, checking these conditions amounts to simply 
an algebraic exercise. A concrete example of this is given in 
the next section. 

The basic idea underlying the analysis in the remainder of 
this subsection is this. Recall that Xjq{k) is defined by ( II 81 . 
Suppose we associate the discrete time k with points on 
the real line spaced apart by a distance proportional to 5n- 
Then, the above technical assumption implies that x^v(fc) is, 
in some sense, close to the solution of a PDE of the form z = 
f(s, z, Vz, V 2 z) with boundary condition h(s, z, Vz, V 2 z) = 
0. Because the Markov chain X/v(fc)/A/ is close to XN{k), as 
established in the last subsection, it is also close to the solution 
of the PDE. The remainder of this subsection is devoted to 
developing this argument rigorously. 

Fix T > 0. Assume that there exists a function z : [0, T] X 
V -)■ R that solves the PDE 

z(t, s) = f(s, z(t, s),Vz(t, s), V 2 z(t, a)), (24) 

with boundary condition 

h(s, z(t, s), Vz(t, s)V 2 z{t, s)) = 

and initial condition z(0, s) = zq(s). Here, V l z(t, s) repre- 
sents all the ith order partial derivatives of z(t, s) with respect 
to s, where i = 1, 2. 
Define 

dt N = S N /M. (25) 

Define 

K N = \T/dt N \ and t N (k) = kdt N . 

Define 

z N (k,n) = z(t N (k),v N (n)) 

and let z N (k) = [z N (k, 1), . . . , z N (k, N)] T G R N . 

Denote the co-norm on R N by || • ||^ • That is, for x G M N , 
with the nth element being x(n), 

IMI^ } = max \x (n) | . 

l<n<N 

Denote the co-norm on R NxK « also by 



That is, for 



x = . . . ,x(K N )] e R NxKn , where for k = 1,...,K N , 

x(k) = [x(k, 1), . . . , x(k, N)] T G R N , we have 

IMl!x^ = ,_ ma x \x(k,n)\. 



k = l,...,K A 



Now we present a lemma on the relationship between the 
z N (k) and f N . 

Lemma 4: Assume that z is continuously differentiable in 
t. Then for each N, there exists UAr(fc) G M. N such that for 
k = 0, . . . , K N - 1, 

z N (k + 1) - z N (k) = -^fN{zN{k)) + dt N u N (k), (26) 



WnW^ = 0(max{7Ar,dtAr}), 



(27) 



where WA r = [u N (0), . . . , u N (K N - 1)] G M ArxK ". 

Proof: Since z is continuously differentiable in t, there 
exists < ci < oo such that for each N, for k = 0, . . . , ifjv — 
1 and n = 1, . . . , N, there exists a function rjy : [0, T] x V — > 
R such that 

z/v(fc + 1j n ) — z^(k, n) 



dt 



N 



z{t N (k),VN(n)) - z(tN(k),v N (n)) 



dt 



N 



z(t N (k), v N (n)) + r N (t N (k),v N (n)), 



(28) 



and for N sufficiently large, rAr(tAr(fc), VN{n))\ < cidt^. 

By (fJTJi and (124b . there exists < C2 < oo, such that for 
each N, for k = 0, . . . , if at — 1 and n = 1, . . . , N, there exists 
a function cf> N : [0, T] x P -> R such that 

z(t N (k),v N (n)) 

= f(v N (n),z N (k, n),Vz N (k, n), V 2 z N {k, n)) 

= fN(z N (k), v N (n))/S N + (f> N (t N (k), v N (n)), (29) 

and for N sufficiently large, \4>N{tN(k), VAr(n))| < C2JN, 
where {7jv} is as defined in (1221 ). 

For each N, for k = 0, . . . , K N - 1 and n = 1, . . . , N, let 

u N (k,n) = (j) N (t N (k),v N (n)) + r N (t N (k), v N (n)), 

and Uiv(fc) = [u]v(fc, 1), . . . , u N (k, N)] T £ R N . Then there 
exists < c < co such that for each N, 



H^iViioo 

Hence ftTJ\ follows. 



By d28]l and 

n = 1, . . . , N, 



^ < cma,x{~/ N ,dt N }. 
for each N, for k = 0, . . . , ifjv — 1 and 



zjv(fc + 1) - z N (k) /jv(zjv(fc)) 



+ u N {k). 



dt]y (5/v 

By this and (ED, we have j26l . ■ 
In the following we show that under some conditions, 

XN(k) and ZAr(fc) are asymptotically close for large N. 
For each N, for k = 0, . . . , A" at and n = 1, . . . , N, define 



£Ar(fc, n) = z^{k, n) — x^{k, n), 



(30) 



and let e N (k) = [e N (k, 1), . . . , e N (k, N)] T G R N . 

By ( fT8l , ( T2o*l ), and ( |30l , we have that for each TV, for A; = 
0, . . . , -ZVa?, there exists UN(k) as defined in Lemma [4] such 
that 



£Ar(fc + 1) = e N {k) + ^(f N (z N {k)) - fN{xN{k))) 



dtNUN{k). 



(31) 



Suppose that for each N, /at G C 1 . Let D/n(x) be the 
derivative matrix of the function fjy at x. Then we have that 
for each N, for k = 1, . . . , Kn and n = 1, . . . , JV, there exists 
a function /jv : R w -> R w such that 



fN(z N (k)) - f N (x N (k)) = Df N (z N (k))e N (k) + f N (e N (k)) 
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and 



Mo) = o. 



(32) 



Then we have from PTT ) 



e N (k + 1) = e N {k) + — (Df N (z N (k))e N (k) 
M 

+ f N {e N (k)))+dt N u N (k). (33) 

Further suppose that for each N, 

11^(0)11^=0. (34) 

Define £Ar = [s N (l), . . . , s N (K N )} G WL N x Kn . Then by 02}, 
(l33l . and (l34l . for each N, there exists a function iJ^ : 

m jvxk„ R w^iv such that 



EAT — Hn(un)- 

It follows that Fat(0) = and £ C 1 . 
For each AT, define 



fXjsr = lim sup 

<*->0 .. (jv) 



Lemma 5: Assume that 



\H N (u)\\W 



(JV) 



z is continuously differentiable in t; 

for each N, f N G C 1 ; 

for each TV, ( 134b holds; and 

the sequence {/iat} is bounded. 



Then 



Ml 



(JV) 



0(\\u N 



Proof: By definition, for each jV, there exists <5 > such 
that for a < 6, 



sup 



|#jv(u) 



II" 



< /i/v + 1. 



0. Then there exists Nq and 



By 03, as N -> oo, 

ax such that for N > jVo, < o>\ < S. Hence, for 

N > N , 

||#jv(u)||W 



< sup 



|#jv(«) 



|(JV) 



OO 



< Mat + 1. 



,l«ll~' \\u\\^>< ai \\u\\\ 

Therefore, there exists < c < oo such that for N > No, 

= l|ffjv(«Jv)||W < (mjv + 1)||^||W 
<(c+l)\\u N \\W. 

This finishes the proof. ■ 
Lemma states that as N — > oo, |[ejv||oo 0, and at least 
with the same rate as ||itjv||oo^. 

Let X N = [X N (l)/M,...,X N (K N )/M], x N = 
[x N (l),...,x N (K N )], and z N = [z N (l), . . . , z N (K N )}, all 
in G M. NxKn . Now we present the main convergence theorem 
of this paper, which states that the value of the normalized 
Markov chain at time k and node n, is close to that of z at 
the corresponding point (t^(k), u/v(n)) G [0, T]xV for large 
M and N. 



Theorem 1: Suppose that the assumptions in Lemma Q] and 
Lemma |5] hold. Then 

II^JV - znW^ = 0(ma,x{'f N ,dt N }) a.s. 

Proof: By (l27t and Lemma [5] there exists < Co < oo 
such that for N sufficiently large, 

ll £ Jv||^ < c max{'j N ,dt N }. (36) 

Let f N in Corollary □ be T/S N . Then K N := [f N M\ = 
[T/dtpf\ := Kn. Hence by Corollary Q] for any sequence 
{Cjv }, for each N, we can take M sufficiently large such that 

oo oo 

J2 p{\\x N - x N \\yp > c/vi < E l i N2 < w - 



JV=1 



JV=1 



(35) By the first Borel-Cantelli Lemma 

p(lim sup {\\X N -x N \\W >C/v}l =0, 

L JV^oc J 

which implies that, a.s., for N sufficiently large, 

\\X N -x N \\W < Cjv . 

Take (at such that for N sufficiently large, 

0v < max{7/v, dt N }. 
Then by the triangle inequality 

\\X N - z N \\W < \\X N - Xtf \\Vp + \\x N ~ z N \\W 
= \\X N -x N \\W + \\e ff \\W, 
a.s., there exists < c < oo such that for N sufficiently large, 

\\Xn - zn\\£' < cmax{'j N ,dt N }. 

This finishes the proof. ■ 
This theorem states that as M — > oo and N — > oo, Xn 
converges uniformly to zn a.s., and at least with the same 
rate as max{'jN,dtN}- 

D. Convergence of Continuous-time-space Extension 

In the following we study the convergence of the 
continuous-time-space extension of the Markov chain X^(k) 
to the PDE solution. Set TV = T/Sjf. For each N, we can 
construct X N(t) and x n(c) with time interval of length 
1/M, with t G [0,Tjv]. Respectively, let X pN (t) and x pN (t), 
where t G [0, T], be the continuous-space extension of X N\t) 
and x Ar(f) (with t G [0, TV]) by piecewise-constant space 
extensions on T) and with time scaled by Sn so that the time- 
interval length is Sn/M := dt^. By piecewise-constant space 
extension of X n, we mean that we construct a piecewise- 
constant function on T) such that the value of this function at 
each point in T> is the value of the component of the vector 
X n corresponding to the grid point that is "closest to the left" 
(taken one component at a time). Then for each t, X p ]y(t) 
and x P N(t) are real-valued functions defined on V. Fig. |2] 
is an illustration of a; at and x p m in a one-dimensional case. 
For fixed T, both X p M{t) and x P N{t) with t G [0, T] are in 
the space D v [0,T] of functions of [0, T] x V -> K and are 



s 



x^kn)&x p Jt,s) 



n & s 




k&t 



• x N {k,n) <^ -^x pN {t,s) 

Fig. 2. An illustration of zjv and x p jv in a one-dimensional case. 

Cadlag with the time component. Define the oo-norm 
on D v [0,T], i.e., for x G D v [0,T], 



i(p) 



|(p) - 



sup |x(i, s)|. 



ts[o,T], 

.ED 



First we show that x p m and z are asymptotically close for 
large N . 

Lemma 6: Suppose that the assumptions in Lemma [5] hold. 
Then 

\\x P N — z||ro — 0(max{7Ar,diAr,c!sAr}). 

Proof: For each TV, for fc = 0, . . . , T^at and n = 1, . . . , TV, 
by the definition of x p n, we have that XpN(tjy(k) 7 viy(n)) = 
XN{k,n). Let f2jv(fc,n) be the subset of [0, T] x 2? con- 
taining (tjv(^)) Wjv(ri)) where x p jv is piecewise constant, i.e., 
(tN(k),VN(n)) £ flN(k,n) and for all (t, s) e QAr(fc, n), 
x P ]sr(t,s) = x P N(tN{k),VN(n)). (For example, for £> C K, 
Uzv(A;, ra) = [ijv(fc), tjv(fc + 1)] x [«jv(n), UAr(n + 1)].) Then 
for each TV, 



\KN-z\\£<\\e N \\W 



max 

k = 0,...,K A 



sup 

(t,s)ef2jv(fc,n) 



|2r(tjv(A;), VN{n)) - z(t,s)\. 



Since z(i, s) is continuously differentiable in t on a compact 
domain, it is Lipschitz continuous in t. Similarly, it is Lipschitz 
continuous in s. Hence there exist < c±,C2 < oo such that 
for each TV, 



max 

= 0,...,Kj 
^=1 > ...,N 



fc=0,...,Kjv, ( i>s ) e n w (fc,n) 



sup \z(t N (k),v N (n)) - z(t,s)\ 

sup \\(t N (k),v N (n)) - (t, s) 

r{k,n) 

< C2 maxjcZs tv , eft at } , 



< c\ max 

^T::^' (t,«)enj,(*,n) 



where || ■ || is some norm on [0, T] x V. Hence, by this and 
there exists < c < oo such that for TV sufficiently large, 

\\ x pN ~ z\\^ < cma,x{-f N ,dt N ,ds N }. 



This finishes the proof. ■ 
Now we present a convergence theorem for the continuous 
functions. 



Theorem 2: Suppose that the assumptions in Lemma Q] and 
Lemma |5] hold. Then 

\\XpN-A\°o =0{m-Ax{-f N ,dt N ,ds N }) a.s. on [0,T] x V. 

Proof: By Lemma [3] , for any sequence {Cat}, for each 
TV, we can take M sufficiently large such that 

oo oo 

£ p{\\x oN - XoN \\g> > c^v} < E 1 / n2 < w - 



jV=l 



N=l 



By the first Borel-Cantelli Lemma 



P<j lim sup {\\X oN -x oN \\£ > ( N } J- = 0, 



AT- 

which implies that, a.s., for TV sufficiently large, 

\\X oN - x oN \\t ] <(won [0,f N ]. 

Since X p m and x p m are the piecewise continuous-space exten- 
sions of X n and x n by constant interpolation, respectively, 
it follows that for any sequence {Ca/}> we can take M 
sufficiently large such that, a.s., for TV sufficiently large, 

\\X p n - XpN\\£ <Cn on [0,71x5. 

Take £jv such that for TV sufficiently large, 

Cat < max{7AT, dtw, dsjy}. 

Then by the triangle inequality 

\\X pN - < \\X pN - x pN \\£ + \\x pN - z||W 

and Lemma [6] a.s., there exists < c < oo such that for TV 
sufficiently large, 

\\X p n — z\\$ < cmax{7jv, dtjy, dsN} on [0, T] x V. 

This finishes the proof. ■ 
This theorem states that as M — > oo and TV — ► oo, the 
continuous-time-space extension X p n of the Markov chain 
XN{k), converges uniformly to z, the solution of the PDE 
a.s., and at least with the same rate as max{7Ar, dtN, dsw}- 
The solution of the PDE can be found quickly by mathemat- 
ical tools readily available and then be used to approximate 
the Markov chain Xi^(k). We give an example of this in the 
next section. 



IV. Application to the Modeling of Large 
Networks 

In this section we present an example of the application 
of our approach to network modeling. We show how the 
Markov chain representing the queue lengths of the nodes in 
the network can be approximated by the solution of a PDE 
using the results of the preceding section. 
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X N (k + 1, n) - X N (k, n) = < 



1 + G(k, n), with probability 

(1 - W(n,X N {k,n)/M)) 
x [P r (n - l)W[n - l,X N (k, n - 1)/M)(1 - W{n + l,X N (k, n + 1)/M)) 
+ Pi(n + l)W(n + l,X N (k, n + 1)/M)(1 - W(n - 1, X N (k, n - 1)/Af ))]; 
— 1 + G(k, n), with probability 
W(n,X N (k,n)/M) 
x [P r (n){l - W(n + 1, X N {k, n + 1)/M))(1 - W(n + 2, X N (k, n + 2)/M)) 
+ P;(n)(l - VK(n - 1, X N {k, n - 1)/Af))(l - VK(n - 2, n - 2)/M))]; 

G(k,n), otherwise. 



(37) 




Fig. 3. An illustration of a wireless sensor network over a two-dimensional 
domain. Destination nodes are located at the far edge. We show the possible 
path of a message originating from a node located in the left-front region. 

A. Network Model 

We consider a network of wireless sensor nodes uniformly 
placed over a domain. In a random fashion, the sensor nodes 
generate data messages that need to be communicated to the 
destination nodes located on the boundary of the domain, 
which represent specialized devices that collect the sensor 
data. The sensor nodes also serve as relays in the routing 
of the messages to the destination nodes. Each sensor node 
has the capacity to store messages and decides to transmit or 
receive messages to or from its immediate neighbors at each 
time instant, but not both. This simplified rule of transmission 
allows for a relatively simple representation. We illustrate 
such a network over a two-dimensional domain in Fig. [3] 
The communication is interference-limited because all nodes 
share the same wireless channel. We assume a simple collision 
protocol: a transmission from a transmitter to a neighboring 
receiver is successful if and only if none of the other neighbors 
of the receiver is a transmitter, as illustrated in Fig. @] 

B. Continuum Model in One Dimension 

We first consider the case of a one-dimensional network, 
where N sensor nodes are uniformly placed over a domain 
Pel and labeled by n — 1, . . . , N. The destination nodes 
are located on the boundary of V, labeled n — and n = 
N + 1. Again let ds^ be the distance between neighboring 
nodes. Let X^(k, n) in (TT~5b be the queue length of node n at 
time k. Let M in ( TTol ) be the maximum queue length of each 
node. 

At each time instant k = 0,1,..., node n decides to be 
a transmitter with probability W(n, Xj^{k, ri)/M). Assume 




Fig. 4. An illustration of the collision protocol: reception at a node fails when 
more than one of its neighbors transmit (regardless of the intended receiver). 



• G 



O 




• • • • 

• • • • • 

l 1 — ~ 1 1 1 — ~. 1 1 

1 n-l n n+l N N+l 

Fig. 5. An illustration of the time evolution of the queues in the one- 
dimensional network model. 

that node n randomly chooses to transmit to the right or the 
left immediate neighbor with probability P r (n) and P/(n), 
respectively. Define G{k) = [G{k, 1), . . . , G(k, N)] T , where 
G(k, n) is the number of messages generated at node n at time 
k. We model G(k,n) by independent Poisson random vari- 
ables with mean g(n). The destination nodes at the boundaries 
of the domain do not have queues; they simply receive any 
message transmitted to it and never itself transmits anything. 
We illustrate the time evolution of the queues in the network 
in Fig. |5] 

The sequence X^(k) defined above forms a Markov chain 
whose evolution is described by ( [Tol l. According to the behav- 
ior of the nodes, the nth component of Fjv(Xjv(fc)/M, U(k)), 
where n = 1 , . . . , N, is defined by d37b at the top of the 
page, where Xiy(k,n) with n < or n > A^ + lare 
defined to be zero. For simplicity, in the following parts we 
set W(n, X?j(k, n)/M) = Xjv(fc,n)/M, which corresponds 
to the transmission rule that a node transmits a message 
with a probability proportional to its queue length. With this 
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. , N, is 



, Xjvj" , the nth component of 



simplification, for x — [x\, . 
Fm{x, U(k)), where n = 1, . 

1 + G(k, n), with probability 

(1 - x n )[P r (n - l)x n -x{l - x n+ i) 
+ Pi(n + l) 

— 1 + G(k, n), with probability 

x n [P r (n)(l - x n+1 )(l - x n+2 ) 
+ Pi(n)(l - x„_i)(l - x n -2)]; 
G(k,n), otherwise, 

where x n with n<0orn>iV + l are defined to be zero 



Define /jv as in ( fT71 ). It follows that for a; = [ 



xi, . 



.,X N 



the nth component of fisr(x), where n = 1, . . . , N, is 

(1 - x n )[P r (n - l)x n _i(l - X n+ i) 

+ Pi(n + l)x n+1 (l - X n -i)] 

- x n [P r (n)(l - x n+ i)(l - x n+2 ) 

+ Pi{n){\ - x„_i)(l - x„_ 2 )] + g(n), (38) 

where x n with n < or n > JV + 1 are defined to be zero. 
Define the deterministic sequence Xjv(fe) as in ( fT8l . 
Set <5jv, defined in Section UlI-CI to be ds^. Let 

dt N = 8 N /M = ds 2 N /M. (39) 

Assume 

Pi (n) = pi (v N (n) ) and P r (n) = p r (v N (n) ) , (40) 

where pi(s) and p r (s) are real- valued functions defined on V. 
As in Section |n] we again assume 

Pi(s) — b(s) + ci(s)dspf andp r (s) = b(s) + c r (s)dsN. (41) 

Let c = ci — c r . Again we call b the diffusion and c 
the convection. In order to guarantee that the number of 
messages entering the system from outside over finite time 
intervals remains finite throughout the limiting process, we 
set g(n) — Mg p (vj^(n))dt^, Assume 6,c;,c r , and g p are in 
C 1 . Then f N e C 1 . 

Let fN(yN,SN) be defined as in Section ITlI-CI Then we 
have the / in (|2TT >: 

f = b( S )^- s {(l-z(s))(l + 3z(s))z s (s)) 

+ 2(1 - z{s))z s (s)b s (s) + z{s)(l - z(s)) 2 b ss (s) 

+ j-(c( s )z( S )(l - z(s)) 2 ) + g p (s). (42) 

Here, recall that, a single subscript s represents first derivative 
and a double subscript ss represents second derivative. 

Based on the behavior of nodes n = 1 and n — N next to 
the destination nodes, we derive the boundary condition for 
the PDE. For example, the node n = 1 receives messages 
only from the right and encounters no interference when 
transmitting to the left. Replacing x n with n < or n > N+l 
by in (1381 . it follows that the 1st component of Jn{x) is 

(1 - x n )Pi(n + l)x n+1 

- x n [Pi(n) + P r (n)(l - x n+1 )(l - x n+2 )] + g(n). (43) 
Similarly, the A^th component of Jn(x) is 
(1 - x n )P r (n - l)x n -x 



- x n [P r (n) + Pi(n)(l - x n _i)(l - x„- 2 )] + g(n). (44) 



Set /3/v, defined in Section BlI-CI to be 1. Then we have the 
h in 



h = -6(a)z(a) 3 + Hs)z(s) 2 - b(s)z(s). 



(45) 



Solving h = for real z, we have the boundary condition 
z(t, s) — 0. This equation might seem confusing to some 
readers as the limit of d43l and (l44i >. if it has not been 
noticed that, unlike /, g is the limit of a different function 

fN(yN,S N )//3 N . 

For fixed T, let z : [0,T] x V -s- R be the solution of 
the PDE ( 1241 . with boundary condition z(t, s) = and initial 
condition z(0, s) = Zo( s )> where the right hand side of 
is 

b(s)-^ ((1 - z(t, s))(l + 3z(t, s))«,(t, *)) 

+ 2(1 - z(t, »))«,(*, s)6 s (s) + 2!(t, s)(l - «(t, s)) 2 fe ss (s) 



+ ^(c( S )z(t, s )(l-z(i, s )) 2 )+ 5p ( s ). 



(46) 



In the following we show the convergence of the Markov 
chain X/v(/c) to the PDE solution z to in the one-dimensional 
network case. Define Km, u n, and en as in Section llH-CI 
Throughout this section we assume 041 holds. By d38l and 
j46l ), it follows that there exists < c < oo such that for N 
sufficiently large, 



hivll^ < cds N . 



(47) 



Albeit arduous, the algebraic manipulation in getting d42l . 
d45l) . and d47T > amounts only to algebraic exercises, the concept 
of which is no more sophisticated than that in getting ( TT4l 
in Section [TT] In practice, we accomplishe such manipulation 
using symbolic tools provided by computer programs such as 
Matlab. 

By (O, for each A, for k = 1, . . . , K N and n = 1, . . . , N, 



we can write £j^(k,n) = H^'' L> (ujy), where 



r(*v 



is 



a real-valued function defined on K ArxKjv . It follows that 



H 



(fc,; 



N 

Define 



(0) = and H N 



(fc,n) 



eC 1 . 



DHn = max 

fc=i,. 

Tl=l 



E 



i=l,... 
J=l.. 



<9iJ 



(fc,n) 
N 



du(i,j) 



(0) 



where is in R NxK ". 

Lemma 7: We have that for each N, 

< DH N . 

Proof: For each N, we have 

r (k,n) 




(0)«(*,J) 



Q H (k,n) 

N tt(0) 



l«(*>i)l 
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= DH N \\u\\W. 
Thus, for each N, for all u^O, 



DH N > 



maxfc=i,...,jsr JV 

n=l,...,N 



Ml 



(48) 



For each N, let w = [f(l), . . . , v(Kn)], where v(k) — 
[v(k, 1), . . . , v(k, N)] T , where for k = 1, . . . , Kn and n = 
1,...,N, 



v(k, n) — sgn 



(feo,no) 



du(k, n) 



-(0), 



where 



{k ,n ) G argmax £ 



k = l,...,K N 

n=i,...,N 't 1 ; -■-'- ff ™ 



(fe.ra) 
N 



du(i,j) 



(0) 



Then 



maxfc=i,...,K JV 



AT 



ll« 



I (AO 



By this and (l48l l we have 



-D-ff/v = sup ■ 



maxk=i,...,K N 

n=l....,N 



I (AT) 



(49) 

By Taylor's theorem, for each N, for k = 1, . . . , Kn and 
n = 1, . . . ,N, there exists H^' n \u) such that 



dH 



N 
N 



j = l,...,N 



(o)tt(y) + 4 M W, (50) 



and for i = 1, . . . , Kn and j = 1, . . . , N, 

&(fc.«V„., 

0. 



ti-s-0 



Hence for each e > 0, there exists <5 such that for 
we have 



<5, 



H { N k ' n) (u)\ 
HI 



(N) 



< e. 



Then for lldl^ < a < 8, 



\H { N k ' n) (u)\ 
sup — - — — — < e. 



| (AT) 



Therefore, for i = 1, . . . , Kn and j = 1, . . . , N, 
lin, sup l^^=0. 



(51) 



By (|50j, for each TV, 

PM«)ll£°< 



max 

h=l,.,.,K N 



H ( N k ' n \u) 



max 

k = l,...,K N 
n = l,...,N 



£ ^ ,n) 5u(<,j)(o)«(*,j) 



<=i * JV 

j = l,...,JV 



Hence 



fJ-N < lim sup 

Halloo s a 



maxfc=i,...,jf N 

n = l r ...,N 



H 



(k,n) 
N 



maxfc=i,...,K N 

ri=l,...,JV 



l«ll 



(W) 



/ 



Hence by d49l ) and ( TSTl ), we finish the proof. ■ 
Notice that DHn is essentially the induced oo-norm of the 
linearized version of the operator Hn- 

Now we present a lemma on the condition of the sequence 
{^n} being bounded for the one-dimensional network case. 

Lemma 8: In the one-dimensional network case, assume 
that the function 



maxi \z , \z. 



\z sa \,\b\,\b s \,\b sa \,\c\,\cs\} (52) 



of (i, s) is bounded on [0,T] x V. Then {/ijv} is bounded. 
Proof: Define 

A N {k) = I N + -±jDf N (z N (k)), 

where In be the identity matrix in l w x N . It follows from 
that for each N and for k = 0, . . . , if at, 

e w (fc + 1) = A N (k)e N (k) + /jv(£ ^ (fc)) + dt N u N {k). 



M 



It follows that 



e N {k) = dt N {A N (k - 1) . . . A(l)u N (0) 
+ A N (k-l)...A(2)u N (l) 



+ 



Define 



D 



(k,n) 
N 



+ A N (k - l)u N (k - 2) + u N (k - 1)) 
+ jj(A N (k-l)...A(2)f N (s N (l)) 
+ A N (k-l)...A(3)f N (s N (2)) 



+ A N {k-l).f N (e N (k-2)) 
+ f N (e N {k-l)). 



0, 0<n<k-3; 
In, n = k — 3; 

A N (k-l)...A N (n + l), n>k-2. 

(53) 



It follows that 



dHp n \u). (fcii ) 



Hence by Lemma [7] 

Mat 



< max V 



•6=1, ...,7V i=l,...,K N 
j=l,...,N 



dt 



N- 



(54) 
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By 438}, for fixed N, for x = [. 
component of Df^{x) 



dx rt 



xi, . . . , xn] t , the (n, m)th 
-(x), where n, m = 1, . . . , N, 



Pi(n)x n (l - ar„_i), 

(1 - a;„)[P r (n - 1)(1 - 
-P/(n + l)x n+1 ] 
+Pl(n)x n (l - x n _ 2 ), 

-[P r (n - l)a; n _i(l - x n+1 ) 
+Pi(n + l)x n+ i(l - x n -!)] 
-[P r (n)(l - x„+i)(l - x n+2 ) 
+P/(n)(l - - ar n _2)], 

(l-x„)[P(n + l)(l-a;„_i) 
-P r (n - l)a: n _i] 
+P r (n)x n (l - x n+2 )}, 

P r (n)x n (l - x n+ i), 





to = n + 1 ; 
rn = n + 2; 
other wise, 



where x n with 7i < or n > N + 1 are defined to be zero. 

Denote the induced oo-norm on R NxAr again by || ■ \\^\ 
That is, for A e R NxN , with the (i, j)th element being 



N 



which is simply the maximum absolute row sum of the matrix. 
Then we have, 



\A N (k) 



\(N) 



1 



™ ax N — (\Pi(n)z N (k,n)(l - z N (k,n- 



l))l 
-1)) 



+ |(1 - ZAr(fc, u))[P r (n - 1)(1 - z N (k,n 

- Pi(n + l)z N (k,n + l)} 

+ Pi(n)z N (k,n)(l -z N (k,n-2))\ 

+ \M- [P r (n- l)z N (k,n-l)(l - z N (k,n + l)) 

+ Pj(n + l)z N {k, n + 1)(1 - z N (k, n - 1))] 

- [P r (n)(l - z;v(fc, n + 1))(1 - ZAr(fc, n + 2)) 
+ P(n)(l - ZAr(fc,n - 1))(1 - ZA r(fc, n - 2))]| 
+ |(1 - z N (k, n))[Pi(n + 1)(1 - zjv(fc, n - 1)) 

- P r {n — l)z7v(fc, Ji — 1)] 

+ P r {n)z N (k 1 n)(l — z N (k,n + 2))]| 
+ |P r (n)zjv(fc, n)(l - ZAr(fc, n + 1))|). 

Put d40j>, (ST}, and the Taylor's expansions (Qj}, (O, and 
(fT~3T > of z,b, and c, respectively, into the above equation and 
rearrange. (Again we omit the detailed algebraic manipulation 
here.) Then we have that there exist < c\ < oo such that 
for each N, for k = 1, ... , Kn and n — 1, . . . , N, 



\\A N (k)\\W 

< max | - c s (v N (n)) - b ss (v N (n)) 

n—l,...,N 

- 2b(v N (n))z ss (t N (k),v N (n)) 
+ 4b ss (v N (n))z(t N (k), v N (n)) 
+ 2b s (v N (n))z s (t N (k),v N (n)) 
+ 4c s (v N (n))z(t N (k),v N (n)) 



+ 4c(v N (n))z s (t N (k),v N (n)) 
+ 6b(v N (n))z s (t N (k),v N (n)) 2 

- 3b ss (v N (n))z(t N (k), v N {n)) 2 

- 3c s (v N (n))z(t N (k) 7 v N (n)) 2 

+ 6b(v N (n))z(t N (k),v N (n))z ss (t N (k),v N (n)) 

- 6c(v N (n))z(t N (k), v N (n))z s (t N (k), v N (n))\ 

ds 3 1 

ds^ ds^ 
:= max \ q (t N (k),v N (n))\— + Ci— + 1. 

n=l,...,N M IVl 



Since d52l is bounded, there exists < c 2 < oo such that 
\q(t,s)\ < c 2 for all (t,s) e [0,T] x P. Hence for each iV 
and for fc = 0, . . . , i£jv, 



ds^ 
~M 



\\A N (k)\\w<i + c 2 ^ + Cl 'r^ 



dsj 



Hence there exists < C3 < 00, for N sufficiently large and 
for k = 0, . . . , Kn, 

dt 2 

\\A N (k)\\W < l + c 3 ^ = l + C 3 dt N . 



Hence by C133D and (154l l. for N sufficiently large, 

UN < K N dt N {\ + c 3 dt N ) KN . 

Since T < 00, there exist < C4 < 00 such that for each TV, 
KjydtN < C4. But as N — > 00, ifjv — > 00, and 



(1 + c 3 di 



-> e 



c 3 T 



Therefore, there exist < C5 < 00 such that for each N, 
fi>N < c 5- This finishes the proof. ■ 
Proposition 1: In the one-dimensional network case, sup- 
pose that the assumption in Lemma [8] holds. Then 

II^JV - Zivll^ = 0{ds N ) a.s. on [0,T] x V. 

Proof: By (|39| > and d47] i. there exists < c < 00 such 
that for N sufficiently large, 

max{-fN , dsN , dtN} < cdsN- 

One can now easily verify that the assumptions in Theorem [TJ 
hold. Then by Theorem [TJ the desired result holds. ■ 
This proposition states that in the one-dimensional network 
case, as M — > 00 and N — ► 00, Xjy converges uniformly to 
zjv a.s., and at least with the same rate as dsN- Analogously, 
for the continuous-time-space extension X p n of X^{k), given 
the same assumption as in the above theorem, by Theorem [2] 
we have 

ll^piV - z\\<£; = 0(ds N ) a.s. on [0,T] x V. 
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z(t ,s) 



N=20, M=N J 



ds N ds N 



d.Sy 



dsn ds N 



dSy 



Fig. 6. The PDE solution z(t, s), at t = t a approximating the normalized 
queue lengths of a one-dimensional network. 



1) Interpretation of the approximation PDE: Now we make 
some remarks on how to use a given approximating PDE. First, 
for fixed N and M, the normalized queue length of node n 
at time k, is approximated by the value of the PDE solution 
z at the corresponding point in [0,T] x V, i.e., 



z((t N {k),v N {n))) R 
Second, we show how to interpret 



X N (k, n) 
M 



C(t ) 



z(t a , s)ds N , 



the area below the curve z(t a ,s) for fixed t a £ [0,T]. Let 
k Q = \t /dtN\. Then we have 



z(t , v N (n))ds N 



X N (k ,n) 
M 



ds 



N, 



the area of the nth rectangle in Fig. [6] Hence 



C(t ) 



N 

£ 

n=l 



N 



z(t ,v N (n))ds N 



n=l 



X N (k ,n) 
M 



ds 



N- 



the sum of all rectangles. If we assume that all messages 
in the queue have roughly the same bits, and think of ds^ 
as the "coverage" of each node, then the area under any 
segment of the curve measures a kind of "data-coverage 
product" of the nodes covered by the segment, in the unit 
of "bit-meter". As N — > oo, the total normalized queue 
length 2~2 n =i XN{k , n) /M of the network does go to infinity; 
however, the coverage ds n of each node goes to 0. Hence the 
sum of the "data-coverage product" can be approximated by 
the finite area C{t ). 

2) Comparison between PDE approximation and Monte 
Carlo simulation: One dimension: We compare the PDE 
approximation obtained from our approach with the Monte 
Carlo simulations for a network over the domain T) = [ — 1,11. 

2 

We use the initial condition Zo(s) = l\e~ s , where l\ > 
is a constant, so that initially the nodes in the middle have 
messages to transmit, while those near the boundaries have 
very few. We set the message generation rate g p (s) — lie~ s , 
where l 2 > is a parameter determining the total load of the 
system. 




-1 -0.5 

o Monte Carlo simulation 



0.5 1 
- PDE solution 



Fig. 7. The Monte Carlo simulations (with different N and M) and the PDE 
solution of a one-dimensional network, with 6=1/2 and c = 0, at t = Is. 



We use three sets of values of N = 20, 50, 80 and M = N 3 , 
and show the PDE solution and the Monte Carlo simulation 
results with different N and M at t = Is. The networks have 
diffusion coefficient 6 = 1/2 and convection coefficient c = 
in Fig. [7] and c = 1 in Fig. [8] respectively, where the x-axis 
denotes the node location and y-axis denotes the normalized 
queue length. 

For the three sets of the values of N = 20, 50, 80 and M = 
N 3 and with c = 0, the maximum absolute errors of the PDE 
approximation are 5.6 x 10~ 3 , 1.3 x 10~ 3 , and 1.1 x 10~ 3 , 
respectively; and with c = 1, the errors are 4.4 x 10~ 3 , 1.5 x 
10 -3 , and 1.1 x 10~ 3 , respectively. As we can see, as N 
and M increase, the resemblance between the Monte Carlo 
simulations and the PDE solution becomes stronger. In the 
case of very large N and M, it is difficult to distinguish the 
results. 

We stress that the PDEs only took fractions of a second 
to solve on a computer, while the Monte Carlo simulations 
took time on the order of tens of hours. We could not do 
Monte Carlo simulations of any larger networks because of 
prohibitively long computation time. 

C. Continuum Model in Two Dimensions 

Generalization of the continuum model to higher dimen- 
sions is straightforward, except for more arduous algebraic 
manipulation. Now we consider the two-dimensional network 
of N\ x N2 sensor nodes. The nodes are uniformly placed 
over the domain T) C K 2 and labeled by (n, m), where 
n = 1, . . . , Ni and m — 1, . . . , iVjj. Again let the distance 
between neighboring nodes be ds^. Assume that the node at 
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N=20, M=N J 




-1 -0.5 

o Monte Carlo simulation 



0.5 1 
PDE solution 



Fig. 8. The Monte Carlo simulations (with different N and M) and the PDE 
solution of a one-dimensional network, with 6 = 1/2 and c = 1, at t = Is. 



location (n, m) randomly chooses to transmit to the north, 
east, south, or west immediate neighbor with probabilities 
P e (n,m) = bi(s) + c e (s)ds N , P w (n,m) = bi(s) + c w (s)dsN, 
P n {n,m) = b 2 (s) + c n (s)ds N , and P s (n,m) = b 2 (s) + 
c s (s)dsN, respectively. Define c\ — c w — c e and c 2 = c s — c n . 

The derivation of the approximating PDE is similar to those 
of the one-dimensional cases, except that we now have to 
consider transmission to and interference from four directions 
instead of two. We present the approximating PDE here 
without the detailed derivation: 



2 



d 



3 = 1 



dz 
ds. 



dbj 



d 



. _^ +z( l_ z) 4^ + - ( (1 _ z) 4) 

dsj ds, ds- 



with boundary condition z(t,s) = and initial condition 
z(0, s) = Zq(s), where t G [0, T] and s — (si, s 2 ) G T>. 

1) Comparison between PDE approximation and Monte 
Carlo simulations: Two dimensions: We compare the PDE 
approximation and the Monte Carlo simulations of a network 
over the domain D = [—1,1] x [—1,1]. We use the initial 
condition Zo(s) = lie~( Sl+S2 ' > , where l\ > is a constant, 
so that initially the nodes in the center have messages to 
transmit, while those near the boundary have very few. We 
set the message generation rate g p (s) — l 2 e~( Sl+s z\ where 
l 2 > is a parameter determining the total load of the system. 

We use three different sets of the values of Ni x N 2 and 
M, where N t = N 2 = 20, 50, 80 and M = Nf. We show 
the contours of the normalized queue length from the PDE 




x 10"' 



x 10" 



-0.5 0.5 
Monte Carlo simulations 



0.5 0.5 
PDE solution 



Fig. 9. The Monte Carlo simulations (from top to bottom, with Ni = 
N 2 = 20,50,80, respectively, and M = Nf) and the PDE solution of 
a two-dimensional network, with bi = 62 = 1/4 and c\ = C2 = 0, at 
t = 0.1s. 



solution and the Monte Carlo simulation results with different 
sets of values of Ni, N 2 , and M, at t = OAs. The networks 
have diffusion coefficients bi — b 2 = 1/4 and convection 
coefficients c\ = c 2 = and ci = —2, c 2 = —4 in Fig. [9] and 
Fig. [TO] respectively. It took 3 days to do the Monte Carlo 
simulation of the network at t — 0.1s with 80 x 80 nodes and 
the maximum queue length M = 80 3 , while the PDE solved 
on the same machines took less than a second. We could not 
do Monte Carlo simulations of any larger networks or greater 
values of t. 

For the three sets of values of N± = N 2 = 20, 50, 80 and 
M = Nf and with ci = c 2 = 0, the maximum absolute errors 
are 3.2 x 10~ 3 , 1.1 x 10~ 3 , and 6.8 x 10~ 4 , respectively; and 
with ci = — 2,c 2 = —4, the errors are 4.1 x 10~ 3 , 1.0 x 
10 -3 , and 6.6 x 10 -4 , respectively. Again the accuracy of the 
continuum model increases with Ni, N 2 , and M, 

V. Conclusion and Future Work 

In this paper we analyze the convergence of a sequence 
of Markov chains to its continuum limit, the solution of a 
PDE, in a two-step procedure. We provide precise sufficient 
conditions for the convergence and the explicit rate of the 
convergence. Based on such convergence we approximate the 
Markov chain modeling a large wireless sensor network by a 
nonlinear diffusion-convection PDE. 

With the sophisticated mathematical tools available for 
PDEs, this approach provides a framework to model and simu- 
late networks with a very large number of components, which 
is practically infeasible for Monte Carlo simulation. Such a 
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Fig. 10. The Monte Carlo simulations (from top to bottom, with N\ = 
N 2 = 20, 50, 80, respectively, and M = N'() and the PDE solution of a 
two-dimensional network, with b\ = = 1/4 and ci = — 2,C2 = —4, at 
t = 0.1s. 



tool enables us to tackle problems such as performance anal- 
ysis and prototyping, resource provisioning, network design, 
network parametric optimization, network control, network 
tomography, and inverse problems, for very large networks. 
For example, we can now use the PDE model to optimize 
some performance metric of a large network by adjusting 
the placement of destination nodes or the routing parameters 
(coefficients in convection terms), with relatively negligible 
computational overhead compared with that of the same task 
done by Monte Carlo simulation. 

The approximation approach can be extended in future 
work with more specific considerations regarding the network, 
which can significantly affect the derivation of the continuum 
model. For example, we can seek to establish continuum 
models for other domains such as the Internet, cellular net- 
works, and traffic networks; we can consider more boundary 
conditions other than sinks, including walls, semi-permeating 
walls, and their composition; the nodes could be nonuniformly 
located, even mobile; transmission could happen between 
nodes that are not immediate neighbors; and the interference 
between nodes could behave differently in the presence of 
power control. 
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